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Abstract 

We show how to construct the optimum superstatistical dynamical model for a 
given experimentally measured time series. For this purpose we generalise the su- 
perstatistics concept and study a Langevin equation with a memory kernel whose 
parameters fluctuate on a large time scale. It is shown how to construct a synthetic 
dynamical model with the same invariant density and correlation function as the 
experimental data. As a main example we apply our method to velocity time series 
measured in high-Reynolds number turbulent Taylor-Couette flow, but the method 
can be applied to many other complex systems in a similar way. 



1 Introduction 



Many complex systems in nature exhibit a rich structure of dynamics, de- 
scribed by a mixture of different stochastic processes on various time scales. It 
has become common to call such dynamical processes with time scale separa- 
tion superstatistical [TR, 6,5,3,2,7,8]. A fascinating feature is that many highly 
nonlinear equations, for example the Navier-Stokes equations, do effectively 
produce stochastic processes that can be very well approximated by linear 
stochastic differential equations with fluctuating parameters on a large time 
scale, i.e., superstatistical processes. Examples where superstatistical tech- 
niques have proved to be a powerful tool are Lagrangian trajectories of tracer 
particles in fully developed turbulence [US], share price dynamics with volatil- 
ity fluctuations [7f8f6] . atmospheric and geophysical processes [9], random ma- 
trix theory [10] and random networks [11]. The basic idea underlying most of 
these models is that locally the system is described by a simple Gaussian 
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stochastic process, and only on a large scale the mixture of various Gaussian 
processes with fluctuating variance leads to more complex behavior. 



So far the superstatistical model building has mainly concentrated on repro- 
ducing the correct stationary density of the nonlinear system under considera- 
tion. However, it is certainly important to incorporate dynamical information 
from the set of all higher-order correlation functions. In fact, experimentalists 
are often able to measure two relevant quantities, the stationary probability 
distribution of a fluctuating time series produced by the complex system un- 
der consideration, and the two-point correlation function for relatively small 
time scale differences. The purpose of this article is to describe how to con- 
struct the optimum superstatistical model for a given complex system taking 
these two experimentally measured quantities as input. As a main example, 
we will construct a suitable superstatistical model for measured velocity time 
series in turbulent Taylor-Couette flow, which is in excellent agreement with 
the experimental data. Our techniques presented here are quite generally ap- 
plicable to large classes of complex systems where the local dynamics is well 
approximated by a generalised Langevin equation with memory kernel. 



2 Standard Langevin equation 

We start from a simple example that will be subsequently generalised. Consider 
locally a linear stochastic differential equation (without memory kernel) of the 
form 



where 7 and (3 are constants and g(t) is Gaussian white noise. This generates 
the Ornstein-Uhlenbeck process, a Gaussian Markov process, with stationary 
distribution p(v) and correlation function c„(A): 



(.) denotes an ensemble average. In the superstatistical approach, one regards 
the parameter j3 as a random variable that fluctuates on a large time scale. If 
this random variable has the probability distribution /(/?), then the stationary 
(marginal) distribution of the complex system under consideration is given by 






C(A) = lim (v(t)v(t + A)) = \e~^ A . 



(2) 




(3) 
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All kinds of distributions P(v) can be generated in this way, for example 
power-law decays in v or stretched exponentials, depending on the properties 
of f{/3) In fact, an experimentally measured P(v) can be used to construct 
f(f3) from eq. (jSJ). The superstatistical correlation function is given by 

C V (A)= d(3f((3)c v (A) = e~^ A / <W(% oc c v (A). (4) 
Jo Jo p 

Clearly, an experimentalist may ask whether such a simple model is appropri- 
ate for a given experimental data set. In most cases it will not, because the 
measured decay of the correlation function is not simply given by exp(— 7 A) 
as is the case for the Ornstein-Uhlenbeck process. In fact, experimental data 
often show local correlation functions that are well described by a sum or dif- 
ference of two exponentials, or sometimes there is oscillating behavior. Our aim 
in the following is thus to extend the local dynamics in such a way that this 
more refined local dynamics, reproducing the correct (measured) correlation 
function, can be incorporated into the superstatistical modelling approach. 



3 Generalised Langevin equation 



A more general starting point is to consider locally a set of coupled linear 
stochastic differential equations 

J^itt) = ~J2 Q j,kVk(t) + E jgj (t). (5) 

The elements of the matrix and the vector H are constants. The noise terms 
are denoted by gj{t). The idea of the theory of Mori-Zwanzig [T2fl~3] is that 
some of the variables Vj(t) are irrelevant, or inaccessible to the experimentalist. 
By eliminating these variables out of the equations of motion, one ends up with 
a set of non-Markovian equations in the relevant variables only. In this article 
we study in detail a 2-dimensional version of this problem, namely 



jv{t) = -jv(t) + V(t) + ag(t) (6) 
jV(t) = - V V(t) + av(t). (7) 

Here, 77 and 7 are positive constants. The constants a and a can have either 
sign. g(t) denotes Gaussian white noise, i.e., (g(t)) = and (g(t)g(t')) = 
8(t — t'). For the above dynamics, the two eigenvalues of the matrix are 
A±B, with 



2A = r/ + 7 and 2B = J(tj - 7) 2 + 4a. (8) 
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It is easy to see that the real parts of these eigenvalues are always positive 
provided 777 — o > 0. In the remaining part of this article, we assume that this 
inequality holds. This insures that the system will locally reach a stationary 
state [12]. 

Our physical interpretation is as follows: The velocity variable v(t) decays with 
relaxation constant 7 and is subjected to two forces: A rapidly fluctuating 
force g(t) that drives the dynamics and another force V(t) that decays with 
relaxation constant rj. This force V(t), however, is not independent of the 
velocity v(t). Rather, the change of V(t) is proportional to the velocity v(t), 
with proportionality constant a. The above linear coupling scheme is one of the 
simplest, but non-trivial, examples of coupled local dynamics between pairs of 
variables. Now assume that an experimentalist can only measure one of these 
two variables. We chose v(t) to be this variable and proceed by eliminating 
V(t) out of the set of equations (|6|7|) . In this way one obtains an equation in 
v(t) with a memory kernel [T21fT5] : 



To obtain this formula, we ignored the contribution V(0) exp(— rjt) because 
this term is unimportant for long times. Expression contains an exponen- 
tial memory kernel and is known in the literature as the generalised Langevin 
equation [T2][T3] . Usually, one assumes that the time dependence of the mem- 
ory kernel and the correlation function of the stochastic term are equal. This 
equality is called the fluctuation-dissipation relation of the second kind [H] 
and the stochastic term g(t) is then interpreted as internal noise. However, 
for the example studied in the present paper, the aforementioned time depen- 
dencies are not equal. Then, the stochastic term g(t) is interpreted as exter- 
nal noise. Such noise terms have also been studied recently in the context of 
stochastic resonance [15] and anomalous diffusion |16j . 

Two well-known special cases are included in our local model. Taking the 
limit r) — » 00 in (Q results in the standard Langevin equation ([1]) without 
memory kernel. Equivalently, this case also arises by putting a = 0, i.e. no 
coupling between velocity and force. Moreover, expressions (jS|7|) reduce to the 
equations of motion of a Brownian harmonic oscillator for 77 = and a < 0. 
The variables v(t) and —V(t) are then interpreted as the velocity and the 
position of the oscillator, respectively. 



4 Local model: stationary distribution 

Before implementing superstatistics, we first have to solve the above local 
model. Let us start with calculating the local stationary distribution p(v) 



d_ 

dt 



v(t) = -jv(t) + a [ dt'e-^-^vif) + ag(t). 



Jo 



(9) 
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Fig. 1. (Colour online) Plot of the correlation function for different values of a. 
The solid line is a representation of expression (|2ip . i.e. the case 5 = 0. The 
dashed-dotted curve shows an oscillatory behavior ([22]) . The dashed and dotted 
curve are respectively a sum and difference of two exponential functions (123j) . 



of the relevant variable v(t). The Fokker-Planck equation [T7j for the joint 
probability distribution p(v, V, t) is 







d d d a 2 d 2 



p(v,V,t). (10) 



The stationary distribution p(v, V) is obtained by setting dp(v,V,t)/dt = 
and is given by 



p(v, V) ~ exp 



7 + V 



a 2 a 2 



[av - r]V f + (777 - a) V 2 



For the marginal distribution p(u) we integrate out the dependence on V 

7 



f +OC O la 2 

p(v)= dVp(v,V) = \ ^-e~ Pv , 

J-00 V 27T 



(12) 



where the parameter is given by 

2(7 + 77) (777- a) 







a 2 (?77 — a + rf) 



(13) 



We see, as required for superstatistics, that the local stationary distribution 
is a Gaussian distribution. 



5 Local model: correlation function 



To proceed further into a dynamical setting, we have to calculate a formal 
solution for v(t). Let us denote the Laplace transform f(w) of a function f(t) 
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by 



dte~ wt f(t) 



From eqs. (ISIITI) one obtains v(w) as 

V(0) 



v[w) 



g(w) 



w + r] 



+ v (0) + ag(w) 



(14) 



(15) 



where the function g(w) is 

q(w) = (w + rj)(w + A + BY l {w + A- 



(16) 



(see ([8]) for the definition of A and B). By calculating the inverse Laplace 
transform of ffl~5D. one obtains 



v{t) = v{0)e{t)+V{0) / dt' e {t-t')e-^' +a dt' g(t - t')g{t'). (17) 



The stationary correlation function c v (A) is obtained by evaluating the en- 
semble average of the product of v(t) and v (t + A) for t — > oo. All the terms 
involving the initial conditions v (0) and 1^(0) that contribute to c v (A) are van- 
ishing for t — y oo and in this limit one ends up with the following expression 
for the correlation function: 



POO 

c v (A) = a 2 / dTg(T)g(r + A). 
Jo 



In order to obtain an explicit expression for c v (A), we have to calculate the 
inverse Laplace transform of g(w) given by ( TL6l) . The result of this transfor- 
mation depends on whether B = or B ^ 0. However it is easy to prove 
that the following inequalities hold: &(A ± B) > given rj > 0, 7 > and 
777 — a > 0. This means that the Laplace transform is well-defined provided 
that the latter inequalities are satisfied. Under the assumption that B = 0, 
one obtains the following expression for the function g(t) 



g(t) = [1 + (rj - A)t]e- At , 
while for the case B 7^ the result is 



(19) 



e (t) = ±(-ri + A + B)e-< A+B > t 



±(- v + A-B)e-^ t . (20) 



We proceed by evaluating (|T8j) to calculate an explicit expression for the cor- 
relation function. Under the assumption that B = 0, one obtains 



-AA 



c v (0) V 2 + A 2 



+ A 2 + A(p 2 - A 2 ) A 



(21) 
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We split the case B^O into two subcases, depending on whether B is imagi- 
nary or real. In case B^O and M(B) = 0, the result for c v (A) is 

^2 = e -AA [T , gin (S / A) + cog (jB / A)] ^ (22) 

C v (v) 

while in case 5^0 and M(B) = B, one ends up with 

c v (A) _ e~ AA 

c v {o) i + r 

with B' = B/i and 



e -BA + Te BA ; (23) 



£' r/ 2 + A 2 + B' 2 ' 5- A(5 + A) 2 -r/ 2 ' 1 ] 

It is easy to show that the signs of T and a are equal. This means that the 
parameter a determines whether one ends up with a sum or a difference of 
exponentials in expression ( 123|) . As a consequence, there are 4 qualitatively 
different types of correlation functions contained in the set of stochastic dif- 
ferential equations (jHE]) • This is illustrated in Figure [TJ The figure shows the 
correlation function for 4 different values of the coupling parameter a but with 
constant values of the other parameters 7 = 0.2 and r] = 0.02. Our chosen 
example values of a are —0.0081 (solid line), —0.1 (dashed-dotted line), 0.001 
(dashed line) and —0.005 (dotted line). Experimentally measured correlation 
functions in turbulence are often of strikingly similar shape, see for example 
Figure 2a in 



6 Superstatistics 



Finally let us apply superstatistics to this local model by introducing large- 
scale fluctuations in (3. We assume that (3, as given by eq. (|T3|) . is distributed 
according to some probability density f((3) and that the changes of (3 occur 
on a time scale that is much larger than the relaxation time of the correlation 
function c v (A). Generally, this results in a superstatistical distribution P(v) 
given by fl3]) that is non-Gaussian. The time-dependence of the superstatistical 
correlation function C V (A) given by @ is equal to the time- dependence of 
c v (A) as long as the /3-fluctuations are produced by fluctuations of a only. 

We are now in a position to formulate the central result of this article. Let 
a complex system be given that exhibits time scale separation and supersta- 
tistical behavior (see [6] for a test). Suppose we know P{v) and C V (A) from 
experimental measurements. How should we construct the optimum synthetic 
superstatistical model that is consistent with the measured data? Of course 
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Fig. 2. (Colour online) Plot of the measured distribution P(v) (•••). The dashed 
line represents a Gaussian approximation 0.15 exp(— 0.07v 2 ) for P(v). The solid line 
is the superstatistical distribution ([3]), with f((3) a lognormal distribution (|25p . 

many nonlinear models are possible, but if one assumes that there is a lo- 
cal linear coupling between pairs of dynamical variables and that locally the 
stochastic process considered exhibits Gaussian behavior, then the model ( joTTj) 
with superstatistical parameter fluctuations is the most appropriate one. The 
distribution f((3) of the parameter (3 can be determined by comparing eq. (J3J) 
with the experimentally measured distribution P(v). Then, depending on the 
shape of the measured correlation function, the values of the parameters 7, 77 
and a can be determined. 

To illustrate our general method, we study as an example time series obtained 
in an experiment performed by Lewis and Swinney [i~9f4] . The data set con- 
tains the values of a single velocity component v(t) as a function of time t 
in turbulent Taylor- Couette flow for Reynolds number Re = 540.000. The 
stationary probability distribution P(v) exhibits non-Gaussian behavior, see 
Figure[2j A suitable choice for f(/3), motivated by the K62 theory of turbulence 
[20] . is a lognormal distribution 

f{fi) = exp (-(ln/3 - hf/2bl) . (25) 

Choosing b\ = 0.42 and 62 = —1.91 one obtains an excellent fit of the ex- 
perimentally measured histogram, see Figure EJ But there are of course many 
different stochastic processes generating the same stationary distribution P(v). 
We can now construct a more precise synthetic model of the dynamics by tak- 
ing into account the information contained in the measured correlation func- 
tion. The experimentally measured correlation function C V (A) / C v (0) is shown 
in Figure [3k, together with a single exponential function exp(— 0.030A) and a 
sum of two exponentials, see expression (123]) with A = 0.024; B = 0.019; T = 
0.571. Figured clearly shows that a single exponential function is not able to 
properly fit the measured correlation function, while a sum of two exponen- 
tials yields an excellent fit. This means that our superstatistical memory-kernel 
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Fig. 3. (Colour online) Plot of the measured normalised correlation function 
C v ( A) jC v (0) (•••). The dashed and solid lines represent an approximation for 
C v (A)/C v (0) in terms of a single exponential exp(— 0.030A) and a sum of two ex- 
ponentials (p3|) . respectively. Fig 3a linear plot, Fig. 3b double logarithmic plot. 



approach is a suitable way to model these data. 

In this paper, we did not consider the /3-process. That is to say, we did not 
construct an explicit model for the time development of the slowly varying 
variable 0. However, it is clear that for very long time scales the details of 
this process will influence the correlation function C v (A)/C v (0) and the sta- 
tionary distribution P(v). To investigate the behaviour for larger time scales 
A, Fig. [3b shows the measured and modelled correlation function in a double- 
logarithmic plot, which emphasizes the long-time behaviour. We see that our 
model (solid line) fits the measured data well up to time scales A w 400. 
Beyond that, there are large statistical fluctuations and the behaviour is dom- 
inated by details of the /3-process. For example, an asymptotic power law decay 
could be constructed by considering a /3-process whose long-term correlations 
decay with a power law. Taking the /3-process into account goes beyond the 
scope of the current paper but is an interesting topic for further research. Our 
current approach as presented in this paper gives testable predictions for small 
and intermediate time scales A. The asymptotic A — > oo behaviour is then 
dominated by the chosen /3-process. 

Let us emphasise that obtaining a sum (and not a difference) of two expo- 
nentials for the correlation function f[2"3"j) is non-trivial. For this it is crucial 
that a > and as a consequence rj ^ and 7 7^ 0. Such correlation functions 
are also observed in biological time series for trajectories of motile cells [21] 
and have been recently obtained theoretically for two-dimensional stochastic 
motion with uncorrelated fluctuations of the speed and the direction of the 
motion |22|. 
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7 Discussion 

To summarise, in this article we have shown how to construct an optimum 
superstatistical dynamical model that exhibits the same invariant density and 
correlation function as an experimentally measured time series extracted from 
a complex system. The relevant universality class of complex systems to which 
our approach is applicable consists of systems where locally pairs of dynam- 
ical variables are coupled in the simplest possible (linear) way, as described 
by eqs. (EE])- This leads to a Langevin equation with memory kernel whose 
parameters fluctuate on a large time scale in a superstatistical way. As an 
example, we applied our approach to experimentally measured velocity time 
series in turbulent Taylor- Couette flow, obtaining excellent agreement with 
the experimental data. Our approach is quite generally applicable to many 
complex systems with time scale separation, as long as high-precision experi- 
mental data on the stationary density and correlation function are available. 
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